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Figure 1: Transitions in the early universe. Going backwards in time there was the decoupling of photons 
(DEC), observed today as the Cosmological Microwave Background, a point of equal energy in 'radiation' 
and (non-relativistic) 'matter' (EQ), a period of nucleosynthesis (NS), probably a QCD transition, and pre- 
sumably an Electro Weak transition. There are good arguments for a period of inflation, with support from 
the microwave background. When inflation ended is uncertain. Leptogenesis (LG) and/or baryogenesis 
(BG) is supposed to have taken place after inflation and before or during the EW transition. 



1. Introduction 

When indicating various transitions that are supposed to have taken place in the early uni- 
verse [jl|] on a time line (Fig. [I]), one is confronted by the fact that what we know with confidence 
goes back only to the time of nucleosynthesis. Was there a QCD transition? Probably, because 
its temperature of 150 MeV is only two orders of magnitude higher than that at the start of nucle- 
osynthesis. The temperature of the conventional electroweak transition is already three orders of 
magnitude higher than that of the QCD transition. Could it have been different? Where should 
one put the end of inflation, baryogenesis? To make progress we need to understand as much as 
possible of the dynamical behavior of a universe, given our knowledge of fundamental interactions. 
This often concerns non-perturbative phenomena. 

This talk is about numerical simulations in field theory, for early times, around and before 
the supposed electroweak transition, using lattice regularization. We will give an impression by 
examples of what has been done, it is not a review. 



2. Inflation, preheating, defects, baryogenesis, thermalization 

Let us take a brief look at various topics in the title of this section. 
Inflation 

Simple models [Q] are formulated in terms of an inflaton field, a, the Friedmann-Lemaitre- 
Robertson- Walker scale factor a, which satisfy the dynamical equations 

.. „ . V 2 a dV(a,---) 

= o + 3Ho ^ + — V — -> 

a L do 



V 1 ( V(y ) 2 „, 
p = T + ^f + V(a,- •) + ••• 

in flat space. Here H is the Hubble rate and p the energy density. When the potential energy 
of the inflaton dominates and its potential is very flat, there is strong damping and 'slow roll', a 
nearly constant Hubble rate and approximately exponential expansion a <=c e Ht . Inflaton models can 
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Figure 2: Small-field (left) and large-field (right) inflaton potentials. 



be divided into two classes: small-field and large-field models (Fig. ^), depending on a rolling 
away or towards the origin, respectively. Other names are 'new inflation' and 'chaotic inflation', 
respectively. The large-field potential shown in Fig. ||] (right) is an example of a so-called hybrid- 
inflation model: 

V{o,(p) = V Q + l -ii 2 a o 2 + g 2 o 2 ^(p-ii 2 ^(p + X{^(p) 2 , (2.1) 

where <p could be a Higgs field. Flatness of the potential requires relatively small jU CT , typically 
jUo- <S /i. During inflation the inflaton field is large, g 2 o 2 — jl 2 > and <p = (read (<p) = 0). Infla- 
tion ends when g 2 o 2 — \l 2 goes negative. Assuming the <p dynamics to be much faster than that of 
the inflaton, we can solve for <p f <p and substitute it back into the potential: <p ' <p — > (/I 2 — g 2 o 2 )/X. 
This gives the trough near the origin shown in the right plot of Fig. ^| (a similar substitution has 
been assumed in the left plot). When the inflaton has gone down the 'waterfall', the universe is 
assumed to be cold 1 and without appreciable baryon- and lepton-number densities, because it has 
just undergone a huge expansion. The energy in the inflaton has to be transferred to other degrees 
of freedom (d.o.f.) in order for baryogenesis to become possible, such that it ends up eventually 
at the standard hot radiation-dominated stage during nucleosynthesis. This heating up (called 're- 
heating') might take too long in a given model, but non-perturbative pre-heating processes may 
have helped to achieve it. 

Preheating 

Rewriting the potential in the form V(o, • • • ) = V(o) + (g 2 G 2 — jJ. 2 )(p^(p H reveals an ef- 
fective (p mass parameter, 

7 7 7 7 

Its time dependence can have drastic non-perturbative effects: 

• (a) oscillates — > <p resonates [Qj; 

• /i 2 ff < — > tachyonic instability & 

Also the a modes can resonate and a is also tachyonic during the time that d 2 V /do 2 < 0. 
These non-perturbative effects can result in large occupation numbers of a and <p, in a limited 
momentum range. Such processes are called pre-heating. 

'This is different in 'warm inflation', see e.g [ffl. 
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Defects 

Preheating processes often result in the formation of defects, depending on the fields involved: 

domain walls, strings, monopoles (see e.g. ||l|]), textures [0], Q-balls, I-balls, oscillons (see 
e.g. 0, ||]), half-knots, Chern-Simons number-densities, . . . 

They may play a role in baryogenesis, thermalization and structure formation. 
Baryogenesis 

The aim is to explain the baryon to photon ratio tig/nyC^ 6.5 10~ 10 from ng = (inflation) and 
B-number, C & CP violation, during particle physics processes out of equilibrium. Many scenarios 
for doing this have been proposed, see e.g. H [l0|]. Here we shall focus on one for which lattice 
simulations have been carried out: Cold Electro Weak Baryogenesis 

i|n|, m 0,111 El- 



In this scenario one assumes that the electroweak transition was tachyonic of the type discussed 



below (|2.1|), and occurred at the end of a period of inflation when the universe was still cold (so the 
electroweak transition was not caused by the falling temperature of the universe.) The electro-weak 
anomaly in the baryon current 

d^B = 3 <?> 1= 7^2 tt/ Vv^ V 

with q the 'topological-charge' density of the SU (2) gauge field, can be integrated over the transi- 
tion to give a baryon asymmetry 

B(t) = ijdt' J d 3 x{q(x/)). (2.2) 

This can be non-zero because the tachyonic transition occurred under the influence of C and CP 
violation. 

Thermalization 

The important question is: how long does it take to thermalize the expanding universe after in- 
flation? The 're-heating' temperature marks the beginning of the radiation-dominated universe. 



Nucleosynthesis constrains this to r r h > 4 MeV Jig]. Thermalization is a non-perturbative problem 
in out-of-equilibrium quantum field theory. 

3. Classical approximation 

Non-perturbative computations of expectation values observables like 

(0(t)) =TrpO{t)=TrpU t (t)OU(t) 

where p is the density operator and U(t) is the evolution operator in real time, are very difficult in 
quantum field theory. The classical approximation can therefore be very useful, since it is possible 
to solve the fully non-linear field equations on a computer (a perturbative analysis is given in 
[17, 18]). To formulate this we can think of a coherent state fll^ ] or Wigner representation [20] 



(0(<M)) = j £>0 c D7r c p c (0 c ,7r c )O(0 c ,7r c 
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where C and n c are classical conjugate variables. If, under suitable conditions, the functional p c is 
positive (generally not true in the quantum regime), then we can use it as a probability for initial 
conditions. The classical approximation then goes as follows: 

1. draw an initial configuration from p c (0 C5 ^c) 

2. solve classical e.o.m. for <p c and n c 

3. evaluate 0(0 C ,7T C ) 

4. average over initial configurations 

One expects a classical approximation to be applicable in cases where large occupation numbers 
dominate, which depends on the observable and on the initial conditions. Let us take a look at 
occupation numbers. For a free scalar field the usual definition is 

n k = (ala k ), a k = — (fifefe + tflk), 
V2(Ok 

where <pk and % are Fourier modes of the canonical variables, and co k is the corresponding energy 
Wk 2 + m 2 ). We can generalize this to time-dependent quasi-particle energies o\ and occupation 
numbers n k ,n k , defined by 

(<M»k> = ( nk+ ^J ^' ^ k7I k) - ( nk+ 27 ^ 7r k)=" k + ^' (3.1) 

Given the correlation functions, we can easily solve for 04, n k and 7%, and see if n k and % are 
larger than 1. 

Consider now a tachyonic transition at t = t c , when ^ s (t c ) = 0. Assume that the couplings 
are weak enough that a gaussian approximation makes sense. The e.o.m. for a component of q> is 
then 



Two cases for which this equation can be solved are Q21| , |22| ] 

A4 = -M 3 (r-? c ), M 3 = -2g 2 o c 6 c , (3.2) 
valid near f c , and the sudden quench [ |23| ] 

Meff = +A t2 > * < <« (3-3) 

= -At 2 , f > tc. (3.4) 

One finds that nk and % grow faster than exponential 2 and n k + 1/2 — n k — > 0, for & < /c max : 



= jix, for the quench. (3.6) 

Quantum (• • • ) dominated by the growing modes can after some time be well reproduced by the 
classical distribution 



exp 



1 ^ / iei 2 \^\ 2 



max 



± V + - 

2 k<k \n k + l/2 + n k n k +\/2-n k 



^ = ^={(O k <K±n k ). (3.7) 



2 For the quench, n k ~ exp \ /I 2 — k 2 (t — t c )]. 
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Figure 3: Inflaton occupation numbers rik at various times. A resonance is visible at t = 876 (t and k in 
VAv-units). From [Q. 



The condition k < fc max can be relaxed provided that the high-momentum modes can be dealt with 
by renormalization. 

So for tachyonic transitions we obtain well motivated initial conditions for a classical approx- 
imation: 

la. at time t\ > t c when 3> 1 and non-linear terms in the e.o.m. are still small 
lb. ('just the half) alternatively, choose t x = t c 

2. choose between & max as in (JO) or (|3.6|), or k max = 'cutoff (all modes) 



3. draw the initial and % from the gaussian ensemble (3.7) 



The alternative lb. is possible because quantum and classical evolution are formally identical for a 



free field. The initial particle numbers are zero in this case, so in (3.7), + l/2±n/t = 1 /2 whence 



the name 'just the half [23]. In case of fc max ^ 'cutoff, the initial modes with k > k max are set to 



zero. Proposals of this kind appeared earlier in the literature [2C] (see also [24]). 

An attractive formulation for numerical simulations is obtained by first discretizing the action 
on a space-time lattice, with spacings a, <C a s , and then deriving the field equations from the 



stationary action principle. This especially useful for gauge fields [25], and it leads to leapfrog- 



type e.o.m. Some milestones in early universe simulations are the computation of the sphaleron rate 



at high temperature pfl, 27], and the unveiling of 're-scattering' effects after parametric resonance 



Q28Q and tachyonic preheating [QJ. In the following some examples will be mentioned of more 



recent work. 

Example: particle numbers in 'New inflation' [ |29| ] 

The inflaton potential is a small field model of the Coleman- Weinberg form 

1 - ,Al,J"\ l \ . 1 ,..4 



Using an equilibrium one-loop potential for non-equilibrium dynamics is questionable, in princi- 
ple, but it may be justifiable by assuming relatively fast dynamics for the non-inflaton d.o.f. The 
parameters are A = 10~ 12 , v = 10~ 3 rap (nip is the Planck mass), typical values for satisfying the 
CMB constraints. The simulation takes the Hubble expansion into account. Initially (<r) « 0, after 
inflation it rolls into minimum of V where it oscillates, the energy decays into a-particles. Figure 
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mt = 23 mt = 24 rat =27 mt =32 




Figure 4: Higgs-field amplitude in x-y slice at various times. From 

|3| (left) shows how inflaton particle-numbers have become very large due to tachyonic preheating, 
and one can see also parametric resonance around time 876. 

By including another field the transfer of energy to other d.o.f. is studied. The potential 
is taken as V = V(o) + \g 2 o 2 § 2 . This is non-tachyonic for - the particle numbers grow 
exponentially because of parametric resonance. Figure || (right) shows the number densities n = 
J for the two species, for two values of the couplings. The resonance is more efficient 

for smaller g , but n$ <C n a and the inflaton does not loose its energy very fast. So in this case 
preheating is not very efficient and the reheating temperature is estimated analytically to be T r h ~ 
y/m P r a ^ H (« 10 7 GeV for g 2 < X). 



Example: mixing in an inflaton-SU(2)-Higgs model 

The potential corresponds to a large-field model 

V = \ii 2 a o 2 + g 2 o 2 ^(p + X{^(p-\v 2 )\ (3.8) 

in which the inflaton is coupled to the Higgs sector of the Standard Model, albeit with somewhat 
unrealistic couplings, X = 0.11/4, g 2 = 2X. The Hubble expansion is negligible for v = 246 GeV 
(H = 0(10 ) GeV). The initial conditions where drawn from a fit to the ^-dependence of the 
distribution at t = t{, with M = 0.18m, m = \[Xv (cf. (3.2)). Figure ^| shows the Higgs field 



in an x-y slice through the lattice at various times after the transition. Intuitively, the non-linear 
scattering of waves causes a good mixing of field modes and shortens the period needed to reach 
the thermalized state. 



Example: equation of state in an inflaton-U(l)-scalar model [ |30| ] 

The equation of state (e.o.s.), i.e. p = p(p), is a basic ingredient in hydrodynamic descriptions 
of the early universe. In simple descriptions the cosmic fluid is often approximated by mixture 
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Figure 5: Parametric plot of the pressures p^, pq and p a vs the energy density p, for times tm 
and two lattice spacings, a s m = 0.5, 1. From 



85-200, 



of various contributions, such as 'radiation' p = p/3, 'matter' p ~ 0, 'vacuum' p = — p, and it is 
interesting to see if and when such an approximation makes sense. The potential in this study is 
the U(l) version of ( |3.8| ). The energy scale is presumed much higher in this case, with parameters 
g 2 = 10~ 2 , X = g 2 /2, jU(j = 4 10 11 GeV, m H = s/2m = 9 10 14 GeV. There are no gauge fields. The 
Hubble expansion is taken into account. Writing 



<P 



re 



the pressures of the inflaton, Higgs and Goldstone bosons are defined as 



Pa 
Pg 
Pa 



lt2_l 

2 ' e 
±r 2 a 2 

2 U 



2\2 



(Vr/a) 2 + \X{r 2 -v 2 ) 
-i(V«/a) 2 , 

; (Vcr/«) 2 + I(^ + gV)cT 2 



Figure ^| shows the emergence of an equation of state and its break-up into various contribu- 
tions. We see how the Goldstone modes approach the radiation e.o.s. p = p/3, whereas the massive 
Higgs and inflaton modes clearly differ from this. The complete pressure is the sum of these com- 
ponents and this study indicates that a definite e.o.s. emerges around tm 85 after the transition. 
Some of the contributions were interpreted in terms of global strings 



Example: strings and hot spots Q31H 

Global strings are an example of (unstable) topological defects which one expects to be pro- 
duced in a tachyonic transition, and which may be relevant for elucidating the mechanism of tachy- 
onic pie-heating. Other possible objects of interest are 'hot spots', places where the inflaton field 
has 'jumped' over the trough in Fig. (right) and ends up again at values \a\ > o c . The model is 
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Figure 6: Left: transparent surface: a = —O c at time tm = 130; strings: |<p| 2 = v 2 /10 at time for which 



spatial average a = for the first time. Right: transparent green surface: \(p\ 
yellow: a = — <7 C ; time tm = 270. From [^lTj. 



\r/l0; blue: a = a c ; 



the same U(l) infiaton-'Higgs' model as in the previous example (with v 2 /2 — ► v 2 ). The parameters 
are g 2 = 10~ 4 , A = 10~ 2 , m = y/lXv = 3 10 15 GeV, and the expansion rate is taken to have the 



constant value H = 3.6 10 3 m. The initialization is at a = o c when y,£s 



-8 2 e 2 - 



m = 0, a la 'just 



the half with k mSLX = 'cutoff. 

Figure |6] shows surfaces where a has fluctuated to ±a c , and tubes of |<p| 2 = v 2 /10, inside of 
which the centers of strings are located at <p = 0. At the later time tm = 270 the strings have mostly 
decayed, but not the 'hot spots'. The authors of remark that the hot spots play an important 
role in the efficiency of tachyonic preheating. 



Example: Tachyonic electroweak quench 

The emergence and development of particle distribution functions in a tachyonic electroweak 
transition has been studied with the SU(2)-Higgs model 

-^SU(2)H = ^ trF ^ + pM < P) t ° M< P " M V<P + X (<P t( P) 2 > 

with fairly realistic couplings, myy = 82 GeV, myi/mw = \/2- The initial conditions for the Higgs 
field were taken to be the quench (3^), 'just the half, with k max = jJ. = /jih/V2 = 0.25/a s . The 



initial gauge field potentials A were set to zero (since they are on the same footing as <Pk with 
k > k max ), and the canonical conjugates A then followed from the Gauss constraint in the temporal 
gauge (Ao = 0) used for solving the e.o.m. The particle distribution functions for the Higgs and 
gauge fields were computed as in ([0|), after fixing the gauge to 'Coulomb' or 'unitary'. 

Figure |7| shows the development of the W -particle numbers. At first the low-momentum modes 
increase exponentially fast due to their coupling to the unstable Higgs modes. At time 7 there is a 
change of behavior in the tail of the distribution. Approximate thermalization was found to set in 
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(3.9) 



Figure 7: Coulomb-gauge W-particle numbers as a function of momentum, time in units of m H l . From [ |32| ] 



from tmn = 30 or so onwards, and not much happened in 50 < tmu < 100. At tmu = 100 there 
is approximate gauge-independence in the transverse W-particle numbers, but the longitudinal W- 
modes settle at a slower rate. Effective Higgs and W temperatures at that time agreed reasonably 
well. Effective chemical potentials jJ. c hem ~ 0.7 (0.9) for W (H) are expected to vanish much later. 



Example: Cold electroweak baryogenesis [15, 33] 



In this simulation the scenario was tested using an effective CP- violating term in the equations 
of motion following from the lagrangian 



SU(2)H ■ 



35 + 1 „ -,,„ 



Here 3 is the number of families and 5 parametrizes effective CP violation, perhaps from physics 
beyond the Standard Model. Since q is a total derivative, q = d^j^, it is customary to rewrite the 



anomaly equation (2.2) in the form 



B(t) = 3{N cs (t)-Ncs(0)), 

with Ncs = / d 3 x j'c S the Chern-Simons number. A useful parameter is furthermore N w , the wind- 
ing number of the Higgs field, since A^cs ~ N w , is expected to hold when fluctuations have suffi- 
ciently diminished. The parameters of the simulation were my/ = 82 GeV, mh/% = 1> v2, v3 
and 2, quench initial conditions with k max = pL. The lattice implementation of CP violation leads to 
implicit e.o.m. which were solved by iteration at a substantially increased time of computation. 

Figure [8] shows an example of the development of (<p^<p), (Ats) and (N w ). We see Higgs length 
at first rising exponentially and then executing a damped oscillation around its finite-temperature 
equilibrium. The CP violation causes an asymmetry in the Chern-Simons number, close to the 
average winding number, which settles earlier. This leads to a baryon asymmetry (4± 
1 ) 10 5 5, which fits observation for the quite reasonable looking value 5 w 1.5 10~ 5 . 
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Figure 9: 3D plot of Chern-Simons number density for times tm\\ = 7, 8, . . . , 15. Red: positive, blue: 
negative. 
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Example: winding defects p4|] 



A more detailed understanding of the production of Chern-Simons number might enable us 
to get more grip on the effects of different forms of CP violation, e.g. by reducing the number 
of variables through modelling. This production appears to be seeded by regions with winding 
number of ps ±1/2, 'half-knots' [34]. Figure |9| shows an example the development of the Chern- 
Simons number density in the tachyonic electroweak quench. Parameters and initial conditions as 
in Fig. [7[ 



Example: kinetic turbulence and scaling §5^, [36|] 

Recently it was demonstrated that concepts of turbulence apply very well to preheating dynam- 



ics [35]. Coefficients of power laws and scaling relations were derived on the basis of Boltzmann 
equations for classical fields, and applied successfully to results of numerical simulations of an 



inflaton field interacting with itself and with another scalar field [35]. In a contribution to this 



meeting, similar behavior was shown to apply to a tachyonic electroweak transition [36]. This 



simulation included the U(l) hypercharge/electromagnetic field, as it was motivated by investigat- 
ing a seeding mechanism for the currently observed magnetic fields on large scales. Parameters: 
my? = 6.15 GeV, m^/myv = 4.65, M(t\ — t c ) = 5 (cf. (p!2|)), gsu{2)/gu{i) as i n Standard Model. 

Figure [R] shows power-law behavior of the variance of the Higgs and inflaton fields, and 
scaling of the W -particle numbers (defined here in terms of the energy spectrum). 




0.01 



Figure 10: Left: turbulent scaling of Higgs ((p(p^) and inflaton (a a) <*t v , v = —2/ (2m - 1 
m = 4 (H). Right: self-similar SU(2) spectrum: n(k,t) = t~^no(kt~P), q = (7 /2)p, p = 1.1/7. 



), m = 3 (<7), 
From I 



4. 4>-derivable approximation 

At larger times the universe is expected to approach local equilibrium if the dynamics is fast 
relative to the Hubble rate. Then the classical approximation breaks down because the energy 
becomes re-distributed over all field modes (equipartition), either leading to lattice artefacts of the 
Rayleigh- Jeans type, or causing the effective temperature to vanish in case of no cutoff. A quantum 
description is needed to be able to describe true thermalization through scattering and damping. 
This has also a non-perturbative features, e.g. a damped oscillation cos (cot) with y = cA 2 
misses essential features when retaining only low orders of the expansion = 1 — cX 2 H . 
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In diagrammatic language, we need to sum an infinite series of diagrams, in particular for two- 
point functions. One possibility is using the hierarchy of Dyson-Schwinger equations. Recently 



so-called <I>-derivable approximations \ pl[ |38[ ] have been used in numerical simulations, and found 
to be able to capture thermalization [|39|]. It is a real-time effective-action method in which the two- 
point function is treated as a basic dynamical variable, in addition to the mean field. The effective 



action is written as [ iO ] 



r(0,G)=5(0)-^TrlnG + ^Tr0M G + o(0 )G)) 

where 5 is the classical action, G is the two-point function, and O is a sum of two-particle irre- 
ducible (2PI) diagrams, 




+ 




+ 




+ 



(4.1) 



with bare vertex functions and dressed two-point functions G. The equations of motion are given 
by 

5H0,G) 8T^,G) 
84>(x) ' SG(x,y) ' 
and after separation into contributions with definite time-ordering, indicated by >, they have a 
non-local form, e.g. for G: 

[d 2 -n 2 -3X^) 2 -3XG(x,x)]G > (x,y) = 2 J dz° J J 3 zIm[£ > (^,z;0,G)]G > (z,j) (4.2) 

-2jT dz°j d 3 zE > (x,z;0,G)Im[G > (z,3;)]. 

Here the selfenergy, 



50 



8G(x,y) ' 



is again given by a series of diagrams, 




+ 




(4.3) 



Interesting aspects of <I>-derivable methods are: 
- they are 'thermodynamically consistent' [37] 
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- global symmetries — ► Ward-Takahashi Identities, 

e.g. translation invariance — > energy-momentum conservation 

- they are renormalizable gl], g|, ^ El' EE EH 

O-derivable approximations are obtained by truncating an expansion of <I>, guided by the number 
of loops, or the order in a coupling- or 1/Af-expansion [38, fffi ]. Putting T on a space-time lattice 



leads to a numerically viable, albeit computationally expensive scheme. The non-locality of fl4.2| ) 
in space can be taken care of (in the usual homogeneous case) by Fourier transformation, which 
gives an e.o.m. for each Fourier mode. The non-locality in time, while in accordance with causality, 
leads to numerically expensive 'memory kernels'. Moreover, these kernels themselves can only be 
computed in simple approximations, e.g. keeping only diagrams of the type given in the first line 



in (4.3); for example, the diagram shown in the second line has two internal vertices, implying two 



four-dimensional summations over space-time, which is too expensive, numerically. 

Initial conditions are given in terms of and G, and the results of solving the e.o.m. can be 
used directly for computing expectation values (O) of observables 0(<f>,G) - there is no need for 
averaging over initial configurations as in classical approximations. 

Example: chiral quark-meson model [ pOp 

The lagrangian of the model is of the Yukawa form 

-jf = fy^d^Y+^diiOd^o + d^'d^n") + \ }x 2 (o 2 + %"% a ) +gy{o + iY 5 T a n a )y, 
where the x a are Pauli matrices, and <I> is truncated to the simple diagram 




(4.4) 



In ref. [ |50| ] the two-point functions are split into real and imaginary parts, G(x,y) ~ F(x,y) — 
ip(x,y)/2, with spectral functions p and 'statistical functions' F, and the discretization is done 
directly for the e.o.m. in terms of components of the fermion two-point function, as in 



1 „PV 



2 Pj Ofiv 
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Figure 12: Left: particle numbers at time t\i = 14; X = 1 /24 (m H = 71 GeV), Lfi = 22.4, a vj ii = 0.7. Right: 
similar for stronger coupling X = 1/4 (m H = 174 GeV), but smaller volume L\i, = 11 .2. From [pip. 



Remarkably, no fermion-doubling was observed. 

Figure |TT](left) shows spatially Fourier transformed Fy(t,p;t,p) vs time for three values of p 
and two different initial conditions. The Yukawa coupling g = 1 . Memory of the initial condition 
gets lost, which can be phrased in terms of a damping time l/")^ amp ~ 15m -1 , with m a 'thermal 
mass'. The approach to the final values takes much longer. Figure [TT] (right) shows the fermion 
particle-distribution function n p at various times in the form ln(— \ + \/n p ), starting from an initial 
condition out of equilibrium. At large times it approaches the line p/T, corresponding to the Fermi- 
Dirac distribution n p = 1/(1 + e p ' T ) for a massless fermion. So this is clear evidence for quantum 
thermalization. Deviations from the FD form due to the interactions appear to be small, even for 
a Yukawa coupling as large as g = 1 . The approach to the FD form can be analyzed in terms of a 



thermalization time 1 /^ 



herm 



95 m- 



Example: tachyonic 'electroweak' quench [ |51| 1 

In this study the scalar sector of the Standard Model is extended to the 0(N) model 

-J8f = \ d^ad^a -\pi 2 <$>a<$>a + \ H<$>a<$>a? , Of = 1, ... ,2V, 
where A °< 1 /N. The functional <1> is approximated to Next-to-Leading-Order in the 1 /N expansion 



[g8|. |49|], after which N is set equal to the SM value 4. 

Figure [lj shows a comparison with the classical approximation for the same model with N = 4 



(without 1/2V expansion). The initial conditions correspond to the quench (p.4|), so vacuum initial 
two-point function G in the quantum <I>-derivable approximation, and 'just the half in the classical 
approximation, with two choices of k max in (3.7): k max = fi ('unstable initialized') and 'cutoff ('all 
initialized'). It is reassuring that at time t\i = 14 the two classical approximations still agree nicely 



with the <I>-derivable one. (Such confirmation has been found earlier in 1+1 dimensions [|52|].) This 
agreement has diminished at time t\i = 100, and it reduces further as time progresses, although 
the general shape of the classical spectra remains similar to the quantal for many hundreds of t\i. 
Taking the <J>-derivable approximation as a benchmark, the k max = 'cutoff implementation of the 
classical approximation appears to be somewhat more accurate in these plots. 
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5. Conclusion and outlook 

Before putting a quantum field out of equilibrium on the computer, approximations had to 
be made. Classical approximations can be well motivated in a limited time regime and they have 
produced many interesting results. The dominating classical modes have a minimum wavelength, 
which makes it feasible to include topological observables such as winding and Chern-Simons 
number densities, and even corresponding terms in the equations of motion, that are notoriously 
difficult to deal with in the quantum theory. We had no time to go into the question of the size 
of topological defects in gauge theories [j53|]. More work is needed on non-abelian gauge fields 



without Higgs fields, in view of the issues raised in [|54]]. 

In the quantum domain <J>-derivable approximations appear to be working well also for late 
times, and the expected features of quantum thermalization have been observed in scalar and 
Yukawa models. The implementation of Nambu-Goldstone bosons is still an issue [55, 56, 48, p6] ]. 
Gauge theories have to my knowledge not been simulated yet with this method. In case of three- 
point interactions (generic in gauge theories and also in scalar field models with non-zero mean 
field), the two-loop selfenergies have contributions that are prohibitively numerically expensive. 
There is the problem of dependence on the gauge-fixing parameter. It has been argued that such 



gauge dependence is acceptable as part of the approximation [57]. For non-abelian gauge theory 
without Higgs fields, the growth of the coupling in the infrared makes it hard to find good trunca- 
tions. 



With fermions there are the usual issues associated with doubling in lattice simulations [|58 



59Q , although this may be avoidable to some extent in non-gauge theories []50|]. Perhaps the need 
to fix the gauge is a blessing in disguise for chiral gauge theories. Real-time simulations with 
fermions might provide an alternative treatment of finite-density. 

Finally, in some cases the Hubble expansion effectively causes the lattice spacing to grow in 
time, which puts limits on the meaningful duration of the simulation. And here is still lots to do in 
numerical simulations of the Big Bang. 
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